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ABSTRACT 

The dynamical consequences of constraining a numerical model with sea surface height data have been 
investigated. The model used for this study is a quasigeostrophic model of the Gulf Stream region. The data 
that have been assimilated are maps of sea surface height obtained as the superposition of sea surface height 
variability deduced from the Geosat altimeter measurements and a mean field constructed from historical 
hydrographic data. The method used for assimilating the data is the nudging technique. Nudging has been 
implemented in such a way as to achieve a high degree of convergence of the surface model fields toward the 
observations. The assimilation of the surface data is thus equivalent to the prescription of a surface pressure 
boundary condition. The authors analyzed the mechanisms of the model adjustment and the characteristics of 
the resultant equilibrium state when the surface data are assimilated. Since the surface data are the superposition 
of a mean component and an eddy component, in order to understand the relative role of these two components 
in determining the characteristics of the final equilibrium state, two different experiments have been considered: 
in the first experiment only the climatological mean field is assimilated, while in the second experiment the 
total surface streamfunction field ( mean plus eddies) has been used. It is shown that the model behavior in the 
presence of the surface data constraint can be conveniently described in terms of baroclinic Fofonoff modes. 
The prescribed mean component of the surface data acts as a “surface topography” in this problem. Its presence 
determines a distortion of the geostrophic contours in the subsurface layers, thus constraining the mean circulation 
in those layers. The intensity of the mean flow is determined by the inflow/ outflow conditions at the open 
boundaries, as well as by eddy forcing and dissipation. 


1. Introduction 

Numerical models of the ocean circulation have un- 
dergone a considerable development in the last 30 
years. In fact, several steps toward a higher degree of 
realism have been undertaken since the first idealized 
numerical studies of westward intensification (Bryan 
1963; Veronis 1966; Holland 1967). These studies were 
performed with barotropic models of coarse horizontal 
resolution in a rectangular domain. A first fundamental 
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improvement in model “realism” has been the in- 
creased horizontal and vertical resolution. The higher 
horizontal resolution has allowed the inclusion of the 
mesoscale eddy field and its interactions with the mean 
flow, processes previously parameterized often in a 
crude fashion. 

A further step in improving the model realism has 
then been the inclusion of a realistic geometry for 
studying the circulation in specific areas of the World 
Ocean. Models of this type can incorporate the effect 
of the actual coastline and bottom topography in the 
area of interest and allow a more straightforward geo- 
graphical correspondence between the model circula- 
tion and the observed ocean circulation. 

The present state of the art in ocean modeling is 
thus represented by eddy-resolving models in a realistic 
domain. Eddy-resolving models describe the evolution 
and statistical equilibrium of a turbulent ocean. There- 
fore, the general circulation of these models includes 
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both characteristics of the mean circulation and eddy 
statistics. 

In the context of this paper, the term eddies is used 
to define deviations of the model fields from the tem- 
poral mean, thus including ringlike structures as well 
as meanders and changes in the stream position. 

Comparisons of some aspects of the model clima- 
tology with observations have shown the sensitivity of 
the model behavior to the geometry, frictional param- 
eterizations, and boundary conditions. These factors 
seem to affect, in an interdependent and complicated 
fashion, the delicate internal processes of eddy-mean 
flow interactions. For this reason models often fail in 
reproducing some basic characteristics of the ocean 
circulation. A typical example of this difficulty is the 
Gulf Stream area, where features like the separation of 
the stream from the coast, its mean path downstream 
to the Grand Banks, the mean eddy kinetic energy level 
and distribution, as well as the southern and northern 
recirculations, are notoriously very difficult to simulate 
accurately. Recent efforts in the modeling of the Gulf 
Stream region (Thompson and Schmitz 1989; Ezer and 
Mellor 1992) have partially succeeded in realistically 
reproducing the stream separation and the eddy kinetic 
energy level. However, they have also shown the sen- 
sitivity of the model circulation to several interdepen- 
dent factors, such as the lateral boundary conditions 
and the inclusion of thermodynamic forcing ( Ezer and 
Mellor 1 992 ) , as well as the presence of a deep western 
boundary current (Thompson and Schmitz 1989) and 
horizontal resolution (Schmitz and Thompson 1993). 
Although models can be “tuned” to be realistic in some 
of their characteristics, a clear understanding of all the 
physical mechanisms in play, as well as the achieve- 
ment of perfectly accurate simulations, are still to be 
considered as work in progress. For the above reasons, 
one can consider the possibility of using data assimi- 
lation techniques to improve the model climatology, 
including characteristics of both the mean field and the 
eddy statistical properties. Data assimilation is a rela- 
tively new topic in oceanography, relative to the long- 
term experience developed in meteorology. In fact, it 
is only with the advent of the ocean satellite missions 
that datasets with a large coverage over synoptic time- 
scales have become available to oceanographers. The 
limitations that ocean models have in their climato- 
logical characteristics as described above suggests the 
use of these datasets not only for prediction purposes, 
as has been the case in meteorology, but also as con- 
straints to the model behavior from a climatological 
point of view. The datasets obtained from satellite mis- 
sions contain information only about the ocean surface. 
Therefore, from the data point of view, numerical 
models can be seen as “dynamical extrapolators” of 
this surface information to areas (e.g., the deep ocean) 
where data are not available. The goal is to achieve, 
through this blending of data and models, a better de- 


scription of the complete, four-dimensional ocean cir- 
culation. 

A considerable amount of work has already been 
done to test different assimilation techniques and to 
determine how effective they are in constraining nu- 
merical models. Most of these investigations have been 
carried out in the context of the so-called twin exper- 
iments, in which the data assimilated are produced by 
the model itself, thus allowing easy verification of the 
degree of success of the data assimilation process. The 
use of real data represents a further step and poses sev- 
eral new issues that have not yet been fully examined. 
In particular, due to the present state of ocean models 
as described above we need to understand the dynam- 
ical implications involved in the process of combining 
data and models that have somewhat different statis- 
tical characteristics in their climatologies. 

In this study we address these issues using a quasi- 
geostrophic (QG) model of the Gulf Stream area. The 
data we assimilate are maps of sea surface height 
(SSH ), which have been obtained as the superposition 
of SSH variability deduced from the Geosat altimeter 
measurements and a mean field constructed from his- 
torical hydrographic data. The surface data have been 
assimilated by using the “nudging” technique, a tech- 
nique implemented in such a way as to achieve a high 
degree of convergence of the model surface field toward 
the observations. 

A preliminary experiment of this type has been pre- 
sented by Capotondi and Holland ( 1993, henceforth 
referred to as CH ). They have shown the effectiveness 
of this surface data constraint in altering the global 
model behavior. A preliminary comparison with some 
available observations seemed to indicate an improve- 
ment in the degree of realism of the model climatology. 
In the present study a deeper analysis of the surface 
data assimilation process is carried out in order to ra- 
tionalize the assimilation results from a dynamical 
viewpoint. Here in Part I we address the following 
question: How does the model respond to the surface 
data constraint; that is, what are the physical processes 
responsible for the model adjustment and for the final 
equilibrium state when surface data are assimilated? 
In particular, since the surface data have been obtained 
as the superposition of a time-mean component and a 
time-dependent component, we want to understand 
the relative contribution of these two components in 
determining the final results. The aim is to obtain a 
physical understanding of the assimilation process so 
that our findings can be generalized to other assimi- 
lation experiments and ways of improving the assim- 
ilation procedure can be devised. In Part II we will try 
to assess the degree of success of the assimilation pro- 
cedure in improving the model climatology by consid- 
ering comparisons with dilferent types of observations. 

The emphasis of the present study is in understand- 
ing the model response to the prescription of a surface 
streamfunction boundary condition or, equivalently. 
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of a surface pressure boundary condition. Therefore, 
in this initial investigation we adopt the approach of 
assuming a complete and uniformly accurate knowl- 
edge of the surface stream function fields. In the lan- 
guage of estimation theory, this study might be defined 
as a study “in control”: To what extent can a surface 
pressure boundary condition control the dynamical 
system under consideration and why? 

The presentation is organized as follows: in section 
2 we describe the model that has been used for the 
experiments, in section 3 we discuss the data, and in 
section 4 we formulate the mathematical framework 
for understanding the model adjustment and final 
equilibrium state when the surface streamfunction is 
constrained to follow a prescribed field. The way a so- 
lution to this problem can be achieved and the char- 
acteristics of such a solution are illustrated with the 
aid of an analytical example in section 5. The assimi- 
lation experiments are described in section 6. Special 
emphasis is devoted to verify the applicability of this 
theoretical framework to the actual problem we are 
considering. To identify the relative effect of the mean 
and eddy components of the surface data in determin- 
ing the characteristics of the model behavior, we present 
in this section two different assimilation experiments: 
in the first one only the surface mean field is assimi- 
lated, while in the second experiment the total surface 
streamfunction field (mean plus eddies) is used. We 
conclude with a discussion in section 7. 

2. The numerical model 

The model is based upon the closed-basin quasigeo- 

strophic model first discussed by Holland ( 1978) and 

now used in many basic studies of eddy-resolved ocean 

circulation. A review of some of those studies is given 

in Holland ( 1 986 ). A description of the physics of that 

model is given in appendix A. The model has five layers 

in the vertical with layer thicknesses of 300, 450, 750, 
1300, and 2200 m, respectively, from top to bottom. 
The values chosen for the reduced gravity at each of 
the four internal interfaces are 1.222 X 10~ 2 , 1.211 

X 10" 2 , 9.448 X 10" 3 , and 4.368 X 10~ 3 m 2 s"\ re- 
spectively, from top to bottom. The QG streamfunc- 
tions for the velocity fields are considered to be at the 
center of each layer, and the transport in each layer is 
given by the streamfunction times its layer thickness. 
The horizontal resolution is Vs degree of latitude and 
longitude, a resolution necessary for resolving the tur- 
bulent processes in the model (Schmitz and Thompson 
1993). The bottom friction coefficient is 10 -8 s -1 , cor- 
responding to a spindown time of approximately 3.2 
years. The coefficient of biharmonic friction is 2 X 10 10 
m 4 * * * * * s _l . This value corresponds to a decay time of about 
22 days at the length scale of the model grid and to 
approximately 1 60 years at length scales of the order 
of 100 km. The model has a realistic coastline, as can 
be seen in Fig. 1 . As is well known, the QG approxi- 
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Fig. 1 . Streamfunction distributions derived from the Bauer-Rob- 
inson climatology for the model layers 1 (a), 2 (b), and 3 (c). Contour 
intervals are 5000 m 2 s _l for layer 1, 4000 m 2 s“‘ for layer 2, and 
1000 m 2 s -1 for layer 3. 


mation is valid only with relatively small amplitude 
variation in bottom topography. The accuracy of the 
model simulation in the case of finite topography, such 
as that existing in some parts of the Gulf Stream area, 
is not known. We have chosen to exclude from this 
initial investigation the uncertainties associated with 
the topographic issue so that the effects of data insertion 
can be isolated from the effects associated with bottom 
topography. A flat bottom, therefore, is chosen. The 
regional nature of the model and its utilization for 
studying an actual piece of the real ocean require that 
the lateral open boundary conditions, as well as surface 
forcing conditions, be carefully considered. As the QG 
model is not very useful for examining local thermo- 
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haline forcing, those aspects of the physical problem 
will not be considered here. For that purpose a PE 
model will be needed. The surface forcing, therefore, 
is given only by the wind stress. Hellerman’s annual 
mean climatological winds have been used (Hellerman 
and Rosenstein 1983). 

The lateral boundary conditions have been con- 
structed so that the inflows and outflows are confined 
to the three upper layers. The mean streamfunction 
fields used for the assimilation experiments have been 
derived from the climatological data of temperature 
and salinity collected by Bauer and Robinson as ex- 
plained in the next section. They are shown in Fig. 

1 . For consistency the transport at the open bound- 
aries in the three upper layers has been derived from 
the same climatological data in all the experiments 
considered in this study. The Bauer-Robinson data- 
set does not resolve the jet entering at the western 
boundary. The latter has been determined by mass 
conservation from the interior geostrophic flow. The 
total transport of the Gulf Stream at the western 
boundary corresponds to the observed value of about 
52 Sv (Sv = 10 6 m 3 s -1 ), which has been measured 
off Cape Fear by Richardson et al. ( 1969). This is 
approximately the location where 80°W intersects 
the North American coast. Thus, the Gulf Stream 
enters the domain as a western boundary current 
somewhat southwest of Cape Flatteras. Another 
problem with the Bauer-Robinson dataset is the 
missing data in the area of the northern recirculation 
gyre, including the model's northern boundary. We 
have chosen to neglect the inflow through the north- 
ern boundary associated with the northern recircu- 
lation gyre, since this could not be determined from 
the data. This choice leads to the development of 
unrealistic features of the mean model circulation in 
all experiments, with and without assimilation, as 
discussed later. 

The inflows and outflows determined from the cli- 
matological data at the northern, eastern, and southern 
boundaries are thus the following: a branch of the 
stream exits the northern boundary east of 43 °W with 
a transport of 7 Sv; the southern boundary intersects 
the subtropical gyre of the climatological data ( see Fig. 
1 ); a transport of 9 Sv enters the domain west of 63° W; 
while 1 5 Sv leave the domain east of the same longi- 
tude. Most of the outflow ( 39 Sv ) occurs at the eastern 
boundary. At this boundary the model streamfunction 
is relaxed toward the climatological values within a 
sponge layer of 1 2 grid points. This allows some degree 
of meandering of the model outflow around the outport 
defined by the climatological values. To prevent re- 
flection of wave energy, a radiation boundary condition 
is also implemented at the eastern boundary. Addi- 
tional details and discussion about the treatment of the 
open boundary conditions are given in Capotondi 
(1993). 


The mean circulation produced by the model when 
the above boundary conditions are used is shown in 
Fig. 2. The three upper layers of the model are initial- 
ized with the climatological fields in Fig. 1 , while layers 
4 and 5 are initially motionless. The temporal mean 
has been defined over a 4-year period. The mean cir- 
culation developed by the model in the three upper 
layers deviates noticeably from the initial conditions. 
In the model solution the Gulf Stream leaves the coast 
north of Cape Hatteras and then flows eastward almost 
zonally. The broad nature of the eastern boundary 
outflow, as well as the lack of a northern recirculation 
gyre inflow, seems to be responsible for the flow ten- 
dency to “fill" the northern half of the domain. This 
is consistent with the “diffuse separation" found by 
Ezer and Mellor ( 1992) when no slope water inflow is 
prescribed at their northeastern boundary. It is also 
consistent with the sensitivity experiments performed 
by Holland and Schmitz (1989, unpublished manu- 
script) with a quasigeostrophic model similar to the 
one used for this study. At the northern boundary a 
“Fofonoff-type" gyre is formed. These features are 
clearly unrealistic and suggest a reexamination of the 
boundary conditions used. However, although im- 
proving ocean models is a necessary, never-ending 
process, this is beyond our purpose here. The aim of 
this study is to analyze how effective assimilation of 
surface data can be in modifying the behavior of a given 
model and why. 

3. The data 

The absolute sea surface topography, which will be 
used to constrain the surface model fields in the assim- 
ilation experiments discussed below, has been con- 
structed by obtaining the sea surface height variability 
from the Geosat altimeter and the mean SSH from 
climatological temperature and salinity data. In fact, 
as is well known, the geoid error in the altimeter mea- 
surements is larger than the oceanic signal that we want 
to measure; so to eliminate this steady error, the tem- 
poral mean is subtracted from the altimeter measure- 
ments. 

The Geosat dataset has been supplied by the ocean- 
ographic group at the Jet Propulsion Laboratory. It 
consisted of 26 months of data, spanning the period 
from 7 November 1986 to the end of December 1988. 
The raw data consist of observations made by the al- 
timeter along the satellite tracks at the times of the 
satellite passage, approximately every 17.05 days. The 
temporal mean has been subtracted from the data and 
the standard corrections (ionospheric and atmospheric 
path delay, inverted barometer, ocean tides, removal 
of blunder points, orbit error corrections, and time- 
mean removal) have been applied. See Zlotnicki et al. 
( 1989) and Holland et al. ( 1991 ) for a discussion of 
these errors. The resulting signals are referred to as 
“sea surface height variability" or “eddy fields" 
throughout this paper. 
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Fig. 2. Mean streamfunction fields in the five model layers 
obtained in the case in which no data are assimilated at the 
surface. The model is forced with Hellerman annual winds 
as well as inflow/outflow at the open boundaries. The time 
average is performed over a 4-year period, (a) Layer 1 : contour 
interval is 10 000 m 2 s~ l . (b) Layer 2: contour interval is 5000 
m 2 s _l . (c) Layer 3: contour interval is 5000 m 2 s _1 . (d) Layer 
4: contour interval is 4000 m 2 s' 1 . (e) Layer 5: contour interval 
is 2500 m 2 s' 1 . 


The approach adopted in this study is to consider 
the model response when complete surface information 
is supplied. Therefore, a statistical interpolation of the 
data onto the model grid at fixed time intervals has 
been performed. The algorithm used for the statistical 
interpolation, which has been carried out in both space 
and time, is the successive corrections method formerly 
used in meteorology (Tripoli and Krishnamurti 1975). 
The description of the algorithm and the considerations 
that led to the parameter choices are given in appen- 
dix B. 

The final interpolated dataset consists of a series of 
SSH maps at 2-day intervals, spanning the period No- 
vember 1986 to May 1988. During the last seven 
months of the Geosat mission (June 1988 to December 
1988), the data coverage degrades considerably. In fact, 
not only are the descending tracks missing in most of 
the domain but the ascending tracks are also very often 


absent, leading to large areas without altimetric infor- 
mation. Therefore, the data after May 1988 have not 
been used. 

The time interval of two days has been chosen be- 
cause we want to assimilate data continuously in time. 
Timescales shorter than a few days do not seem to be 
relevant in the model dynamics (Holland and Mala- 
notte-Rizzoli 1989 ) so that we can interpolate linearly 
between our bidaily maps to obtain data to assimilate 
at every time step (one hour) of our model run. The 
characteristics of the eddy field that emerge from this 
interpolated dataset are discussed in Part II. 

The mean streamfunction field in the upper layer, 
the initial conditions for the first three layers, and the 
inflows and outflows at the open boundaries have been 
computed from the climatological data of temperature 
and salinity analyzed by Bauer and Robinson and pre- 
sented in version VIII of their atlas (Bauer and Rob- 
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inson 1985). The procedure is as follows. First a dy- 
namic height computation is performed from the ref- 
erence level of 1125 m, which corresponds to the 
middle of layer 3 of our numerical model. The inte- 
gration has been carried out up to the depths of 525 
and 150 m, which corresponds to the middle of layers 
2 and 1, respectively. The resulting dynamic height 
fields can be immediately translated into the stream- 
function fields for the first two layers relative to the 
third layer. The reference flow in layer 3 has been ob- 
tained by assuming that its horizontal structure is the 
same as the one obtained for layer 2 and by requiring 
that the total transport of the inflow at the western 
boundary coincides with the observed value of 52 Sv 
(Richardson et al. 1969). The final result for the cli- 
matological streamfunction fields in the upper three 
layers of the model is shown in Fig. 1 . These fields will 
be used as initial conditions for the corresponding lay- 
ers in all experiments performed in this study. Only 
the field for the first layer, however, will be used in 
conjunction with Geosat data to obtain maps of total 
streamfunction, which will be assimilated into the 
model. The model subsurface layers will be free to 
evolve from the initial conditions and eventually de- 
velop a different mean circulation as a result of the 
interactions with the assimilated eddy field. The use of 
initial conditions obtained from a dynamic computa- 
tion was adopted in order to build into the model a 
“realistic” baroclinic structure that could reduce the 
model adjustment time during the assimilation exper- 
iment. The use of climatological data for determining 
the missing mean component of the Geosat measure- 
ments might be called into question due to the different 
duration of the two datasets. This choice should be 
considered as a “working hypothesis” whose conse- 
quences on the assimilation results will be carefully 
analyzed and discussed. However, as shown in CH, the 
characteristics of the upper-layer mean field ( Fig. I a ) 
show a remarkable correspondence with the distribu- 
tion of eddy kinetic energy obtained from the inter- 
polated altimeter data ( see Fig. 5a in Part II ). Features 
of the mean path, like the large curve around the Grand 
Banks and the subsequent splitting into two separate 
branches, are clearly suggested by the eddy field itself. 
Thus, this choice for the surface mean streamfunction 
field appears consistent and sensible. 

4. A theoretical framework 

The surface data have been assimilated by using the 
“nudging” method. This method, formerly introduced 
in meteorology (Anthes 1974), has been used in 
oceanography in several studies of assimilation of sur- 
face data both in QG and in primitive equation models 
(Verron and Holland 1989; Holland and Malanotte- 
Rizzoli 1989; Holland et al. 1991; Malanotte-Rizzoli 
and Young 1992). A detailed review of the method is 
given in Ghil and Malanotte-Rizzoli ( 1991 ). 


The implementation of nudging in the QG model 
used is straightforward. The upper-layer model equa- 
tion is altered by adding a relaxation term in the form: 

dVVt - , 

— - = ‘physics’ - /?(V V> - VVobs). (4.1) 

Ol 

The equations for the lower layers are left unchanged. 
Here “physics” includes the rate of change of vortex 
stretching, the advection of potential vorticity by the 
surface flow, the steady wind forcing, and the bihar- 
monic friction term. The relaxation coefficient R is, in 
general, a function of space and time to account for 
irregular data and data error distributions in space and 
time (Ghil and Malanotte-Rizzoli 1991; Malanotte- 
Rizzoli and Young 1992). Our emphasis in this study 
is to analyze the model response to the prescription of 
a complete and uniformly accurate surface informa- 
tion. Surface data, in fact, are available at each grid 
point and at time intervals frequent enough to allow a 
continuous assimilation in time. Therefore, the nudg- 
ing coefficient can be chosen to be a constant. The 
value we have chosen is R = (0.5 day )"' , corresponding 
to a relaxation time scale of 1/2 day. This value is close 
to the upper limit dictated by numerical stability con- 
siderations and thus constitutes a “strong” nudging. 
In fact, it can be shown by scale analysis (Capotondi 
1 993 ) that the time scale of 1/2 day is at least one order 
of magnitude shorter than the timescales associated 
with the terms in physics. As a consequence of this 
implementation of nudging, the model upper-layer 
streamfunction will become very close to the “ob- 
served” streamfunction \f/ obs . How do the model sub- 
surface fields respond to this surface constraint? Since 
the surface data include a time-mean component and 
an eddy component and we are interested in separating 
the effect of the time-averaged flow from the effect of 
the eddies in controlling the model dynamics, we for- 
mally split all the model variables, the streamfunction 
\p k and the potential vorticity < 7 *, in the same fashion: 

tk = 'i'k + Vk (4.2a) 

Qk = Qk + Qk\ (4.2b) 

subscripts indicate the layer. The temporal mean is 
supposed to be computed over a time interval much 
longer than the eddy timescales. The expression for the 
QG potential vorticity is given by Eq. ( A4) in appendix 
A. In the presence of a strong nudging the time-aver- 
aged version of the model equations (A5) can be re- 
written in the form: 

« fc*. (4.3a) 

5 + = k = 2, 5. (4.3b) 

at 

The upper-layer equation expresses the data constraint 
we are imposing. The surface wind forcing has thus 
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been replaced by a surface pressure boundary condi- 
tion. The equations for the subsurface layers describe 
the slow variation of the mean potential vorticity in 
the presence of mean advection, dissipation by bihar- 
monic friction, and eddy advection of eddy potential 
vorticity. Biharmonic friction is negligible at the scales 
typical of the mean circulation (Holland 1978). The 
analysis of several numerical QG simulations (Rhines 
and Holland 1979; Holland and Rhines 1980; Marshall 
1984; Lozier and Riser 1990) have shown the existence 
of different flow regimes; in some regions the eddy flux 
divergence is large and important in the mean balance, 
while in other regions eddy fluxes are small in the mean. 
In particular, Marshall’s results with a barotropic model 
show that in some areas close to the western boundary 
eddy fluxes are the driving agent of the mean flow 
across q contours. However, in spite of the local im- 
portance of the eddy flux divergence term, his results 
also show that q remains a strong constraint on \p in 
the inertial recirculation area. The mean flow is so 
strong in this area that even a large eddy forcing causes 
only a relatively small deflection of \j/ across q. The 
importance of mean flow advection in shaping the 
mean potential vorticity distribution in the subsurface 
layers is also evident in the numerical model we are 
using for the present study when no data are assimilated 
(Capotondi 1993). In fact, the model domain mainly 
covers a region of intense mean flow. Since we are 
interested here in the large-scale structure of the cir- 
culation, it appears sensible to adopt an inertial ap- 
proach in our theoretical considerations and neglect at 
first order the eddy flux divergence term in Eq. (4.3b). 
We will verify the validity of this inertial assumption 
in section 6, where the results of the assimilation ex- 
periments are qualitatively compared with the predic- 
tions derived from this dynamical framework. At steady 
state the system (4.3) thus becomes 

h ~ £obs (4.4a) 

J(fa,q k )~ 0, * = 2,5. (4.4b) 

The set of equations (4.4b) defines a generalization of 
the Fofonoff problem (Fofonoff 1954) to a baroclinic, 
four-layer “ocean,” with a prescribed “surface topog- 
raphy” given by ^ obs and inflow/ outflow conditions at 
the boundaries. Marshall and Nurser ( 1986) showed 
how to construct analytical solutions to the baroclinic 
generalization of a Fofonoff problem in an idealized 
rectangular ocean bounded by solid walls. In the con- 
text of our experiments, the particular solution of Eqs. 
(4.4b) that satisfies the given lateral and surface 
boundary conditions cannot be determined analytically 
due to the irregular model geometry and to the speci- 
fication of a “surface topography” and boundary con- 
ditions that are not analytically defined. Therefore, it 
is instructive to consider a simple analytical example 
to illustrate how an inertial solution to the problem 
(4.4) can be achieved. In particular, we are interested 


in identifying the influence of the “surface topography” 
^ obs in determining the characteristics of the global so- 
lution. This is considered in the next section. 

5. An analytical example 

The mathematical formalism adopted for this anal- 
ysis follows the approach introduced by Rhines and 
Young (1982; henceforth referred to as RY) in the 
context of the wind-driven circulation. In the case of 
the wind-driven circulation, the wind stress curl and, 
therefore, the barotropic streamfunction are known. 
In the present case the given quantity is the upper- 
layer streamfunction. As discussed in RY the contri- 
bution of the relative vorticity to the QG potential vor- 
ticity is negligible when considering the large-scale 
steady circulation away from lateral boundaries. 
Therefore, we simplify the expression of the potential 
vorticity [Eq. (A4) in appendix A] by neglecting the 





x-“L x=L 


(b) 


(c) 


Fig. 3. (a) Surface streamfunction for the analytical example; 
yp obt describes an anticyclonic flow inside the disk of radius R . (b) 
Geometry of the geostrophic contours in the second layer. The dashed 
line indicates the disk inside which the surface flow is confined, (c) 
Geometry of the geostrophic contours in layer 3. All the lateral 
boundaries are assumed to be closed. The dashed line encloses the 
disk of radius R ; the dotted line defines the area of closed geostrophic 
contours in layer 2. 
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relative vorticity V 2 ^ A . We consider here the simple 
case in which \p obs is given by the anticyclonic flow 
shown in Fig. 3a and described by the expression 

R 2 -x 2 -y 2 ), r<R 

R (5.1) 

0, r>R. 

The surfa ce veloci ty field is confined inside the disk r 
^ R ( r = Vjc 2 T y 2 ) with intensity increasing from zero 
at the center to the maximum value 2\p 0 /R at the pe- 
riphery of the disk. The presence of this surface flow 
distorts the interface between layer 1 and layer 2, pro- 
ducing a circular depression in the interface. We con- 
sider first a two-layer model and then extend our anal- 
ysis to a three-layer model. A continuously stratified 
case for this particular problem has been discussed in 
Capotondi ( 1993) following the considerations of RY 
in the context of the wind-driven circulation. 



a. Two-layer model 

We consider a rectangular domain as shown in Fig. 
3. The model has a flat bottom, specified boundary 
conditions, and prescribed surface flow. Neglecting the 
relative vorticity and considering the fact that the Ja- 
cobian of \p 2 with itself vanishes, the equation for the 
second layer can be written, from (4.4), in the form: 

J(fa, Py + F^fa*) ~ 0. ( 5 . 2 ) 

The quantity we define 

Qi = Py + F 2] fa bs ( 5 . 3 ) 

is a known function. Therefore, the problem (5.2) is 
a linear problem in which the q 2 contours define the 
mean streamlines. A general solution of (5.2) will be 
of the form 


^2 “ ^2(^2)- 

For \p obs given in (5.1 ), the function q 2 is 
^ 21^0 


(5.4) 


Qi = 


R 2 

Py, 


(R + yoi “ - (y “ T 02 ) ), r < R 


The quantity y 0 2 is defined as 

1 &R 2 

} ; 02 = r ' 


2 F 21^0 


r> R. 
(5.5) 


(5.6) 


where \U mdx \ = 2 fa/R is the maximum surface ve- 
locity. Therefore, in order to have closed contours, the 
surface velocity must exceed the phase speed of the 
long baroclinic Rossby waves supported by the present 
model, c p/l , where | c ph | = p/F 2 1 . In fact, in the system 
we are considering, where only one layer can evolve 
freely, these waves are the only agent that can propagate 
the information about the boundary values into the 
interior. However, if the flow is sufficiently strong to 
oppose the Rossby wave propagation, areas isolated 
from the boundaries can be created. The condition 
( 5.7 ) also coincides with the requirement that the basic- 
state potential vorticity gradient changes sign in the 
lower layer (in the case U is a zonal flow) (Pedlosky 
1979) and can thus be interpreted as a necessary con- 
dition for instability. 

Contours of q 2 are shown in Fig. 3b for the case y 02 
- \ R, which occurs when = PR/ F 2 \ • The value of 

fa on the open contours is determined by the boundary 
values. Therefore, if no inflow or outflow is prescribed 
at the boundaries, motion is possible only inside the 
closed contours, as discussed in RY. Inside closed con- 
tours the flow is undefined at this order. This is a case 
of steady-state resonance in which any flow is in prin- 
ciple possible. A weak forcing can produce an arbi- 
trarily large response, which can be limited only by 
considerations of frictional effects. The selection of a 
particular solution in these areas will be obtained, 
therefore, as a balance between forcing and dissipation. 
In the absence of any external forcing, as in our case, 
the only source of energy for the motion within the 
closed contours is given by the eddy flux divergence 
term. As shown by Holland and Rhines ( 1980), the 
component of this flux that is more effective in driving 
the deep gyres is the flux of interface height, corre- 
sponding to vertical propagation of horizontal mo- 
mentum. This justifies the choice, adopted in RY, of 
parameterizing this flux as a downgradient flux of in- 
terface height displacement 

J(fa, Qi) - — V • [kVF 2 \(\^\ ~ fa)], (5.8) 

where k is the diffusion coefficient, generally a function 
of position. If for convenience k is chosen to be a con- 
stant and dissipation is given by bottom friction in the 
form 


Dl = - t V 2 £ 2 , (5.9) 

the relationship between fa and q 2 is found (see RY 
for the detailed derivation) to be linear inside closed 
contours: 


The q 2 contours are straight lines outside the disk of 
radius R and arcs of circle inside the disk. As in the 
problem discussed in RY, closed contours can be found 
if T 02 < R- This condition is satisfied if 

\U mu \>4~. ( 5 - 7 ) 

F 2X 


dq 2 kF 2 i + £* 

We should notice here that the intensity of the flow in 
layer 2 depends on the forcing and dissipation param- 
eters k and and it is an increasing function of k only 
if € is different from zero. If no explicit dissipation were 
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present, dA 2 /dq 2 would attain its maximum value, 
which only depends upon the model density structure: 


dAj. 1 
dq 2 F 2l ' 


(5.11) 


The same result would hold also with a spatially varying 
diffusion coefficient k. 

The general solution for \p 2 over the whole domain 
can be written in the form: 


*2=^2ft + C 2 . (5.12) 

Suppose that the boundary conditions prescribed in 
the second layer are a uniform eastward flow both at 
the western and the eastern boundaries: 

\p 2 = ~U 2 y at x - — L, L. 

In this case we have 


F 2 \ 


A 


g'nH 2 ' 


If F 2 1 is very small, either because the second layer is 
very deep or because the density difference between 
the two layers is very large, the q 2 contours cannot 
diverge noticeably from the 0y contours. Therefore, 
the flow in the second layer can be expected to be 
mainly controlled by the boundary values. From what 
was shown before, the area of closed contours can be 
expected to be associated with the most intense com- 
ponent of the flow. In the following section we thus 
consider the case of closed lateral boundaries and con- 
centrate on the analysis of the vertical structure of the 
area of closed contours. 


b. Three-layer model 


open contours (r > R) 

closed contours (r < R) y 

where r = Sx 1 + (y - y 02 ) 2 and R = R - y 02 . 

(5.13) 

In (5.12) C 2 is a constant chosen so that \p 2 is contin- 
uous at the edge of the closed contours: 



C 2 


0, 




Qi < $R 

(5.14) 

Qi > 0R- 


Different flow regimes with different relationships be- 
tween streamfunction and potential vorticity are found 
within the closed contours and in the region of open 
contours. 

The potential vorticity q 2 is given by 
Qi — Qi F 2 1^2 

= 42(1 - A 2 F 2i ) - F 2l C 2 . (5.15) 

Inside the closed contours, if bottom friction is absent 
(c = 0), q 2 is a constant given by 


From (4.4b) and (A4) the governing equations for 
layers 2 and 3 at the leading order are 

J(fa,0y + / r 21 &b. + F^h) ~ 0 (5.17a) 

J(h* Py + Fnh) ~ o. (5.17b) 


In the second layer the streamline distribution is now 
determined not only by the prescribed surface topog- 
raphy but also by the topography of the lower interface, 
which is a function of the flow in layer 3. However, 
the intensity of the flow can be expected to decrease 
with depth so that at first order the displacement of 
the lower interface in layer 2 can be considered much 
smaller than the displacement of the upper interface. 
In these conditions an approximate solution in the sec- 
ond layer can be obtained as before, and it is given by 
(5.12). In the case of closed boundaries and in the 
absence of any dissipation in layer 2, the constants A 2 
and C 2 are given by 


r o. 


I F 2 \ + F a ’ 


open contours (r> R) 
closed contours (r < R) 


(5.18) 


and 


q 2 = 0/^1+ U 2 (5.16) 

The constant value of q 2 coincides with the value of 
potential vorticity at the northern rim of the gyre as 
determined by the planetary term / QR and by the sloping 
of the interface associated with the uniform eastward 
flow U 2y which is prescribed by the boundary values. 

From (5.3) we can see that the extent to which 
^obs constrains the pattern of the flow in the lower layer 
is strongly dependent on the density structure of our 
two-layer system as expressed by 


C 2 = -0RA 2 . (5.19) 

The corresponding potential vorticity distribution is 
constant inside the closed contours and coincides with 
the value of the planetary component at the northern 
rim of the gyre, q 2 = /3R. We should recall, however, 
that in the model used for the assimilation experiments 
we have biharmonic friction as a dissipative term. 

Consider now the flow that can be expected in layer 
3. From (5.17b) we have 

^3 = A 3 ( q 3 ), 


(5.20) 
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where 

<73 = f3y + F n A 2 q 2 + F 32 C 2 

= 0y+Ft<**-F J2 0RA 2 , (5.21) 

where 

0 = (1 + F 32 A 2 )0 (5.22) 

and 

F = A 2 Fy 2 F 2 \\ (5.23) 

A 2 is given in (5.18). The q 3 contours, as well as the 
\ p 3 contours, can thus be related to the surface topog- 
raphy i/^. The influence of the surface flow in shap- 
ing the mean circulation in layer 3 is measured by 
F = A 2 F 32 F 2 \ . Outside the closed q 2 contours, A 2 is 
zero, so that no surface information can be felt in layer 
3. Inside the closed q 2 contours, F is given by 


p _ Ft, 2 F 7 1 
F 2 1 + F 23 

If II 2 = H 3 = //, we have that 


F 2 3 = F 32 


fl 

g'liH ' 


(5.24) 


Therefore, (5.24 ) can be rewritten in the form 


1 - _L _L 

F F 2 i + F 23 


(5.25) 


so that the magnitude of F is smaller than either F 2] 
or F 2 3 , and it is given by 


F = 


fl 

Hg y 


(5.26) 


where g = g\ 2 T g 2 3 is the reduced gravity associated 
with the total density difference between layer 3 and 
layer 1. If, more generally, // 3 is different from // 2 , F 
can be written 




/§ i 

gnHs j + P2 ~ Pi 


P3 ” P2 


(5.27) 


We can distinguish two limiting cases. 

1 ) (p 3 - p 2 ) $> (p 2 - p!). In this case F ~ F 32 so that 
the penetration of the surface information is only deter- 
mined by the “rigidity” of layer 3. If layer 3 is very deep 
and/or its density is much larger than the density of the 
layer above, the influence of \p obs in determiningthe dis- 
tribution of the q 2 contours can be expected to be neg- 
ligible with respect to the planetary term. 

2) (p 3 - p 2 ) (p 2 - p,). In this case F 

~~ (/o/# 2 i# 3 )- The influence of the surface infor- 
mation still depends upon the thickness of layer 3, but 
now it depends upon the largest density difference, 
which is the one between layer 1 and layer 2. 


From this simple analysis we see that the penetration 
of the mean surface information is tightly linked to the 
stratification characteristics of the area under consid- 
eration. 

Similar analyses can be extended to a continuous 
stratified case. For brevity we do not consider this case 
here. The interested reader can find the derivation in 
Capotondi ( 1993). As in the case of the wind-driven 
circulation considered by RY, the area of closed geo- 
strophic contours becomes smaller and displaced pole- 
ward with depth. In three dimensions the region of 
closed contours has a “bowl-like” shape whose maxi- 
mum depth, corresponding to the depth of influence 
of the surface streamfunction field, is a function of the 
vertical stratification. 

c. Discussion 

From this analytical example we can make the fol- 
lowing points. 

(i) Due to nonlinear effects, the prescription of a 
surface flow ^ obs can constrain the flow structure in 
the subsurface layers. The “depth of influence” of this 
surface constraint is very strongly affected by the ver- 
tical density profile. 

(ii) The intensity of the flow is dictated by the 
boundary values on the open contours and by the 
characteristics of forcing and dissipation processes in- 
side the closed contours. 

(iii) If closed contours are present in a given layer 
as a consequence of the prescription of the surface field, 
different flow regimes can be found: the functional re- 
lationship between streamfunction and potential vor- 
ticity is dictated by the boundary values in the area of 
the open contours and by the characteristics of forcing 
and dissipation inside closed contours. 

(iv) If the eddy flux of interface height is the only 
nonconservative mechanism present, the intensity of 
the flow inside closed contours will not depend on the 
intensity of the eddy field but only upon the model 
density structure. The corresponding potential vorticity 
fields are constant. However, if some explicit dissipa- 
tion is present in each layer, as in the assimilation ex- 
periments that we are going to discuss, the amplitude 
of the solution will depend on the relative strength of 
forcing and dissipation. 

( v) If no-flow conditions are specified at the bound- 
aries, motion is possible only inside closed contours as 
in RY. The driving agent for this flow is the eddy flux 
divergence term. 

We now return to the specific problem under con- 
sideration and try to verify to what extent the above 
theoretical ideas can explain the results of our assim- 
ilation experiments. 

6. The assimilation experiments 

We consider here two different experiments: in the 
first one only the time-mean component of the surface 
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data is assimilated, while in the second experiment the 
total upper-layer stream function is used in the relax- 
ation term. The comparison between the results of these 
two experiments will help identify the relative effects 
of the mean and eddy components in determining the 
characteristics of the flow fields that develop in the 
model subsurface layers. In the first experiment the 
nudging term has been added to the equation for the 
first layer in the form 

- V 2 ^ obs ), (6.1) 

where \j/ obs is the climatological field in Fig. la. The 
relaxation coefficient is R - (0.5 day) -1 . At each time 
step \[/\ is relaxed toward the steady field \[/ obs with a 
very short relaxation timescale. Therefore, any time- 
dependent motion that the model might try to develop 
will be strongly damped. We will see later, in fact, that 
the time-dependent motion found in this experiment 
is extremely weak. From the considerations developed 
in the previous section, eddies can be expected to be 
the forcing agent for the flow inside the closed geo- 
strophic contours in the subsurface layers. In particular, 
they are the only source of vortieity for the motion in 
layers 4 and 5, where no inflows or outflows are spec- 
ified at the boundaries. Therefore, if the eddy field has 
vanishing intensity, only a very weak time-averaged 
flow can be expected in the two deepest layers. 

The numerical simulation has been carried out for 
20 years to allow all the transient processes associated 
with the initial conditions to decay. The initial con- 
dition in the three upper layers is given by the clima- 
tological fields in Fig. 1, while no flow is assumed ini- 
tially in the two deepest layers. We have monitored 
the time evolution of the total kinetic energy in each 
of the five layers to ensure that statistical steady state 
(in this case almost coincident with an absolute steady 
state) has been reached. The “climatology” of this nu- 
merical experiment has been computed by averaging 
the model fields over the last four years of the simu- 
lation. 

In the second experiment the upper-layer stream- 
function ipi is relaxed toward the total “observed” 
streamfunction i p obs , where 

= ^obs + *obs(0. (6.2) 

As before i£ obs is the climatological field in Fig. la, and 
^obs( 0 is the sequence of eddy maps constructed from 
the Geosat data. As described in section 2, the total 
duration of the \[/ obs dataset is 570 days. The initial 
conditions are the same as for the previous experiment. 
Since we are now imposing a time-dependent con- 
straint at the surface, we need to define sensible criteria 
for assessing when the model has adjusted to the ob- 
servations. The evolution of the total kinetic energy 
during the first 570 days of the experiment shows that 
in each of the five layers the energy increases from the 
initial value to a “steady” level during the first 10-20 
days of the simulation. After this short transient the 


level of total kinetic energy remains practically constant 
in each layer. The potential vortieity distributions on 
the other hand show a continuous evolution. Starting 
from the initial distributions associated with our choice 
of initial conditions, potential vortieity is redistributed 
by advection processes. These processes now include 
not only mean flow advection as in the previous ex- 
periment but also eddy advection. While the mean flow 
tries to establish constant values of potential vortieity 
along streamlines, the turbulent eddy field acts as an 
efficient mixing agent, which tends to smooth out the 
q gradients. The evolution of the potential vortieity 
fields is the result of the competition between these two 
processes, and a statistical steady state will be reached 
when mean flow advection balances the eddy mixing 
effect. In Figs. 4a and 4b we show instantaneous po- 
tential vortieity maps in layers 2 and 3, respectively. 
Notice the convoluted distribution of the q contours, 
which is due to eddy advection and is a manifestation 
of the enstrophy cascade. 

Eddy mixing can be expected to depend only on the 
statistical characteristics of the eddy field and not on 
the details of its instantaneous realizations. To allow 
these processes to evolve until statistical equilibrium 
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Fig. 4. Instantaneous potential vortieity fields in layer 2 (top) and 
3 (bottom) at day 96 from the beginning of the assimilation experi- 
ment. In this experiment the total surface streamfunction field 
= lAob* + ^obs is assimilated. Contour interval is 4 X 10"* s" 1 in layer 
2 and 3 X 10 -6 s _1 in layer 3. 
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is reached, we have extended this experiment beyond 
the 570-day duration of our dataset by assimilating the 
same data in a sequence of runs, each of which is started 
from the final fields of the previous one. The total ex- 
periment consists of 20 of these assimilation segments, 
totaling 1 1 400 days or about 3 1 .6 years of spinup time. 
The convergence of the potential vorticity fields toward 
an equilibrium distribution is illustrated in Fig. 5, where 
we show meridional profiles of mean potential vorticity 
in the four subsurface layers. The different curves in 



LATITUDE (DEGREES) 



LATITUDE (DEGREES) 


each diagram refer to the averages over each of the 
570-day assimilation segments: “A” corresponding to 
the first segment, “B’ 1 to the second, and so on. These 
profiles have been computed as averages over 10° of 
longitude centered at 55°W. In layers 4 and 5 the me- 
ridional profiles are dominated by the planetary vor- 
ticity gradients, and only relatively small adjustments 
take place from segment to segment. In layers 2 and 3 
we can see the convergence of these different profiles 
to a meridional distribution showing a large plateau in 



LATITUDE (DEGREES) 



LATITUDE (DEGREES) 


Fig. 5. Meridional profiles of mean potential vorticity at 55°W in the four model subsurface layers during the experiment in which the 
total surface stream function is assimilated. The sequence of surface eddy stream function fields #>&,(/), covering a period of 570 days, is 
repeated periodically 20 times in order to allow the adjustment of the subsurface potential vorticity fields. The curves in each panel refer to 
the averages over each of the 570-day assimilation segments: A corresponding to the first segment and B to the second and so on. 
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its central part, where the potential vorticity has been 
homogenized due to the very effective eddy mixing. 
We can also notice how the value of q tends to decrease, 
from segment to segment, at the northern end of these 
meridional profiles due to advection of low potential 
vorticity anomalies by the mean flow. 

After the 31.6 years of spinup time, variations can 
still be observed from segment to segment in the mean 
streamfunction fields. The rms differences between the 
stream functions corresponding to the last two segments 
are only a few percent in the upper three layers, but 
they can be as large as 40% in layers 4 and 5, where 
they are mainly associated with slight changes in the 
position of the gyres present in the deep mean fields. 
The “climatology” of this model has been computed 
by considering a time average over the last segment. 

We are now interested in comparing the model cli- 
matologies obtained in these two experiments and in- 
terpreting them in the context of the dynamical frame- 
work developed in sections 4 and 5. To that end we 
first try to infer, on the basis of the dynamical consid- 
erations derived from the analytical example, the 
structure of the flow that can be expected in layers 2 
and 3 when the climatological field i^ obs is imposed at 
the surface. We then analyze in detail the mean cir- 
culation and the mean potential vorticity fields in these 
two scenarios: the first one in which the model time 
variability is very weak, and the second one in which 
the model eddy field is much more energetic and con- 
strained to follow the surface eddy observations. 


a. Streamfunction distributions predicted in layers 2 
and 3 


In the first experiment only a vanishingly weak flow 
is found in the two bottom layers, as we will see below. 
This experiment can thus be conveniently defined as 
one with three moving layers over a motionless deep 
ocean, as in our analytical three-layer case (section 5b). 
Therefore, we concentrate our attention on layers 2 
and 3. In the following, the climatological fields for 
layers 2 and 3, corresponding to Fig. lb and Fig. lc, 
will be referred to as yp 2obs and ^ 3obs , respectively. At 
the initial time the ratio between the standard devia- 
tions of ^ 3obs and f obs is 


VZ ('Piobs)ej 
Vz (fobs )Tj 

ij 


13%, 


(6.3) 


where the indices / and j vary over the entire model 
domain. Therefore, at this time, the function 


Qi = Py + F 2l f obs (6.4) 

defines, to a good approximation, the \f/ 2 distribution. 
Contours of q 2 are shown in Fig. 6a. As in the analytical 
example discussed in the previous section, some of the 
q 2 contours in Fig. 6a are closed and do not reach the 




80°W 70°W 60°W 50 D W 40°vf 


Fig. 6. (a) Contours of the function q 2 = &y + F 2 i^obs for the 
assimilation experiments; \f/ ob5 is the surface climatological stream- 
function field that is assimilated. The q 2 contours supply, approxi- 
mately, the pattern of the flow field in layer 2 during the assimilation 
experiments, (b) Contours of the function = fty + Fjj^obs* where 
^ 2 ot« is the climatological streamfunction field computed for layer 2 
from the Bauer-Robinson data. This field is considered here as a 
good approximation for the time-averaged streamfunction distribu- 
tion in layer 2. The <? 3 contours represent the approximate stream- 
function distribution in layer 3. 


boundaries, while some others are connected with the 
boundaries. In the area of the northern recirculation 
gyre, where ^ obs = 0, the q 2 contours coincide with the 
fiy contours and are given by zonal lines. In this area 
these contours suggest a westward flow emanating from 
the eastward flowing jet. In fact, as the surface topog- 
raphy is flat in this area, fluid particles must move along 
latitude circles to conserve their potential vorticity. A 
no-flow condition must be satisfied at the coastline, so 
that some higher-order physics is required there to close 
the circulation. Therefore, we might anticipate the for- 
mation of a boundary jet along the coastline where 
relative vorticity, or friction, will no longer be negli- 
gible. Away from this region, the circulation suggested 
in layer 2 by the q 2 _ contours does not show evident 
discrepancies with ^ 2o bs (Fig- lb). Therefore, we use 
^ 2 obs in (5.17b) to obtain an approximation for the 
structure of the flow field in layer 3. At initial time 
layer 4 is at rest. In the experiment in which only 
f obs is assimilated, the two bottom layers can be con- 
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sidered motionless. Therefore, the lower interface does 
not enter into play in determining the circulation pat- 
tern in layer 3. The latter can thus be described, within 
the approximations made, by the function q 3 given by 

<73 = ft}' + ^32^2obs- (6.5) 

The contours of q 3 are shown in Fig. 6b. Also in this 
case the absence of a northern recirculation gyre flow 
in i/^obs leads to zonal contours of g 3 , which imply a 
westward flow in this area of layer 3. Again, some 
higher-order physics must enter into play close to the 
solid boundary to satisfy the no-flow condition. 

If we compare the flow pattern in layer 3, as given 
by the <? 3 contours, with the flow pattern in layer 2, we 
can notice how the recirculation gyre in layer 3 appears 
tighter and more elongated in the NE-SW direction 
with respect to the recirculation in layer 2, whose outer 
streamlines extend beyond the model southern bound- 
ary. In Fig. 6b the contours south of about 30°N join 
both the eastern and western boundaries. At these lat- 
itudes the stream function values specified at the west- 
ern boundary are constant, and they have only minor 
variations along the eastern boundary. Therefore, only 
a weak flow can be expected south of 30°N. This vari- 
ation with depth of the shape of the recirculation area 
represents the model analog of the results of the ana- 
lytical example, where the area of closed geostrophic 
contours becomes smaller and displaced northward 
with depth. 

We now compare and contrast the mean stream- 
function distributions that are obtained in layers 2 and 
3 during the assimilation experiments with the “theo- 
retical” flow patterns in Fig. 6. 

b. The mean streamfunction fields 

The time-averaged circulation obtained in the first 
experiment is shown in Fig. 7 in all the five model 
layers. Inthe first layer the flow field is essentially the 
same as V'obs* As anticipated, the nudging term repre- 
sents the dominant contribution in the equation for 
the first layer, so that the upper-layer streamfunction 
becomes almost identical to i^ obs . A closer comparison 
between the two surface fields shows that the major 
differences occur in the area of the jet separation from 
the coast, northeast of Cape Hatteras and in the area 
of the Grand Banks, where one of the branches of the 
stream turnsnorthward. At both locations the stream- 
lines in the \p x field tend to “open” toward the coast. 
This deviation of with respect to \p obs is produced 
by the circulation that develops in the subsurface layers 
north of the stream, as we will see in a moment. The 
maximum differences between \[/ x and ^ obs in both 
places are of the order of 1 5%. Consider now the cir- 
culation in the second layer. If the dynamical frame- 
work developed in sections 3 and 4 captures the essen- 
tial physics of this model simulation, the flow pattern 
obtained in layer 2 should be well described by the 


distribution of the q 2 contours in Fig. 6a. Comparison 
between Figs. 7b and 6a shows, indeed, striking simi- 
larities. The structure of the subtropical recirculation 
gyre obtained in this numerical simulation is rendered 
in great detail by the q 2 contours. The three anticyclonic 
cells, which can be observed in Fig. 6a around 70°W, 
58°W, and 42°W, do appear as features of the time- 
averaged circulation in layer 2. The same is true for 
the cyclonic cell predicted by the q 2 contours inside 
the curve of the stream around the Grand Banks. Also 
as anticipated, a westward flow can be observed in the 
area north of the stream. As the stream emerges from 
Cape Hatteras and flows eastward as a free jet, fluid 
particles detach from the stream and move westward. 
As expected, a thin jet forms along the coastline to 
close the circulation. This coastal jet, whose intensity 
tends to increase with latitude, develops instabilities. 
As it tries to follow the irregular and sinuous coastline, 
meanders develop and ringlike structures are shed, 
which remain trapped between the jet and the bound- 
ary. The presence of this coastal jet, not predicted by 
the simplified derivation of q 2 , has the effect of some- 
how distorting the whole flow field so that a precise 
agreement between \[/ 2 and q 2 cannot be found. The 
presence of the flow in layer 3, which has not been 
considered in the derivation of q 2 , is an additional rea- 
son for discrepancies. However, the basic characteristics 
of the flow in the second layer are predicted by the 
structure of q 2y supporting the hypothesis of the inertial 
nature of the circulation. 

Similar considerations can be applied to layer 3 
(Fig. 7c). In this case the streamline distribution 
should be compared with the distribution of the q 3 
contours in Fig. 6b. Also in this case the shape of 
the subtropical recirculation gyre, the presence of 
smaller-scale anticyclonic cells, the development of 
a westward flow in the subpolar area, and the con- 
sequent formation of the coastal jet are features pre- 
dicted by the <? 3 contours. The most energetic com- 
ponent of the circulation is found north of 30°N, a 
latitude which defines the southern border of the re- 
circulation at this depth. 

In Figs. 7d and 7e we show for completeness the 
time-averaged streamfunction fields in layers 4 and 
5, respectively. In both layers, as expected, the flow 
is quite small almost everywhere. The only noticeable 
component of the circulation is found in the prox- 
imity of the northern boundary, where eddies pro- 
duced by instabilities of the boundary jet have rel- 
atively larger amplitudes. Maps of eddy kinetic en- 
ergy show, in fact, values lower than a few 
centimeters squared per second squared in most of 
the domain. Only close to the northern boundary 
isolated values as high as 100 cm 2 s~ 2 can be found, 
clearly associated with instabilities of the jet. This 
appears to be the only area where eddies can drive 
any flow in the two deepest layers. 
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Fig. 7. Time-averaged stream function field obtained in the 
experiment in which only the mean field ^ obs is assimilated 
at the surface, (a) Layer 1: contour interval is 5000 m 2 s _1 . 
(b) Layer 2: contour interval is 2000 m 2 s“ 1 . (c) Layer 3: con- 
tour interval is 2000 m 2 s' *. (d) Layer 4: contour interval is 
2000 m 2 s~ l . (e) Layer 5: contour interval is 2000 m 2 s -1 . 


The “climatological” streamfunction fields obtained 
in the experiment in which both mean and eddy com- 
ponents are assimilated is shown in Fig. 8. We can 
immediately notice the similarity of the flow patterns 
in the three upper layers with the results of the previous 
experiment ( Fig. 7 ) . The surface layer is strongly con- 
strained by the nudjging procedure so that cannot 
deviate much from ^ obs . However, also the circulation 
in layers 2 and 3, although not directly constrained, 
has basically the same structure as in the experiment 
where only \p obs was assimilated. In both layers the 
shape of the recirculation gyres with all their smaller- 
scale features has remained very similar. Some differ- 
ences can be observed in the area north of the stream. 
In Figs. 8b and 8c we can still notice the tendency for 


fluid particles to leave the jet and move westward, but 
now no well-defined boundary jet develops along the 
coastline. The presence of the eddy field is now able 
to supply a potential vorticity input, which allows the 
mean flow to move northward also in the interior, 
without the need of invoking the higher-order physics 
of a boundary layer. However, a northern recirculation 
gyre is still absent. Eddies alone seem to be unable to 
drive this component of the circulation, thus supporting 
some recent findings of the thermodynamic mainte- 
nance of this gyre (Ezer and Mellor 1992). Due to the 
presence of an energetic eddy field, motion is now pos- 
sible also in layers 4 and 5. The most energetic com- 
ponent of the flow is found in the western half of the 
domain and represents the deep expression of the in- 
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Fig. 8. Time-averaged streamfunction fields from the ex- 
periment in which the total surface streamfunction \p obs is as- 
similated at the surface, (a) Layer 1: contour interval is 5000 
m 2 s'*, (b) Layer 2: contour interval is 2000 m 2 s~‘. (c) Layer 
3: contour interval is 2000 m 2 s (d) Layer 4: contour interval 
is 2000 m 2 s ’, (e) Layer 5: contour interval is 2000 m 2 s _l . 


ertial recirculation for this numerical experiment. The 
streamfunction fields in these two layers are rich in 
small-scale features. Notice, in particular, the tendency 
for the formation of elongated zonal gyres. The pres- 
ence of these gyres could be an artifact of the neglect 
of bottom topography. In fact, in the absence of any 
topographic steering, “free” flow tends to develop along 
fiy contours. We will see in Part II that evidence of 
zonal jets has indeed been found in observations of the 
deep flow in this area. 

Although the structure of the circulation in layers 
2 and 3 has not been noticeably affected by the as- 
similation of the eddy field, the intensity of the flow 
has been altered. To illustrate the differences in the 
flow strength, differences that are induced when sur- 


face eddies are assimilated, we show in Fig. 9 merid- 
ional profiles of mean zonal velocity along 55°W in 
all the five model layers. The profiles have been av- 
eraged over 10° of longitude. In each figure the 
dashed line represents the zonal velocity obtained in 
the experiment in which only i p obs is assimilated, 
while the solid line corresponds to the case in which 
the total i/'obs is imposed at the surface. In layer 1 the 
two profiles are almost coincident as a consequence 
of the strong nudging toward the same climatological 
flow field. In layers 4 and 5 the flow is practically 
zero in the absence of eddies. In layers 2 and 3 a 
strong correlation can be noticed between the solid 
and dashed velocity profiles. The position of the 
eastward jet and of the westward return flow has not 
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been significantly altered by the assimilation of the 
eddy field, but the amplitude of the peak zonal ve- 
locities, both eastward and westward, has been en- 
hanced. 

c. Mean potential vorticity fields 

The mean potential vorticity fields obtained in the 
first experiment are shown in Fig. 10. In layer 1 the 
potential vorticity distribution has remained very sim- 
ilar to the initial distribution (not shown). The major 




LATITUOE ( DEGREES 1 


Fig. 9. Meridional profiles of mean zonal velocity at 55°W 
in the five model layers. The solid line corresponds to the 
experiment in which the total surface streamfunction is as- 
similated; the dashed line corresponds to the case in which 
only the mean field is assimilated. 


differences occur in the “subpolar area,” where the 
stretching effect produced by the flow in the second 
layer determines a distortion of the j 8y contours present 
in the initial field. In layers 2 and 3 the potential vor- 
ticity distributions show the advective control of q. This 
can clearly be seen in Figs. 1 1 and 1 2, where the 
streamfunction distributions in layers 2 and 3, respec- 
tively, are directly compared with the potential vorticity 
distributions in the corresponding layers. In both layers 
the q contours reproduce, in fact, the shape of the re- 
circulation gyre with large areas of reduced gradients 
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Fig. 1 1 . Verification of the degree of agreement between stream- 
function and potential vorticity contours in layer 2. (a) Stream function 
contours. Contour interval is 3000 m 2 s'. (b) Potential vorticity 
contours: contour interval is 2.5 X 10 -6 s '. 


the differences introduced in their potential vorticity 
distributions by the presence of the surface eddy field 
are of particular interest. The maps in Figs. 13b and 
13c no longer show any closed contours reproducing 
the exact shape of the recirculation gyres as in Figs. 
10b and 10c. Only tongues of low potential vorticity 
anomalies can be noticed on the southern flanks of the 
recirculation gyres. As before, the mean flow tries to 
create uniform distributions of potential vorticity along 
streamlines, but now the energetic eddy mixing pre- 
vents the completion of this process. In Figs. 14 and 
15 we compare, as before, the ip contours with the q 
contours in layers 2 and 3, respectively. Although the 
potential vorticity distributions show clearly the effect 
of mean flow advection, the almost exact correspon- 
dence observed in Figs. 1 1 and 1 2 is now lost. 

7. Discussion and conclusions 

In this paper we have tried to understand the mod- 
ifications induced in the model fields by the assimila- 
tion of surface data whose climatological characteristics 
are different from the climatology of the unconstrained 
model. In particular, we have analyzed the relative ef- 
fect on the model behavior of the two components of 


the surface observations: the mean component and the 
eddy component. If the relaxation time in the nudging 
term is much shorter than the typical model timescales, 
the surface fields become practically coincident with 
the observations and can be considered as given. In 
these conditions we have shown that the structure of 
the subsurface circulation is mainly determined by the 
characteristics of the surface mean field. This structure 
is the result of the model nonlinearities and can be 
interpreted in the framework of “baroclinic Fofonoff 
modes” in a domain with a prescribed surface topog- 
raphy (the surface mean field, which is assimilated at 
the surface) and inflow /outflow conditions specified 
at the open boundaries. We have shown with the aid 
of an analytical example how such a solution can be 
achieved. The geometry of the geostrophic contours in 
each layer can be related to the pattern of the clima- 
tological field imposed at the surface. Some of these 
contours join the boundaries, while others are closed 
and isolated from the boundaries. On the open con- 
tours the intensity of the flow is determined by the 
boundary values. Inside the closed contours, on the 
other hand, the amplitude of the flow can be expected 
to be the result of a balance between forcing and dis- 
sipation. Therefore, different flow regimes may exist. 




Fig. 12. Verification of the degree of agreement between stream- 
function and potential vorticity contours in layer 3. (a) Streamfunction 
contours. Contour interval is 3000 m 2 s" 1 . (b) Potential vorticity 
contours: contour interval is 2.5 X 10 -6 s _l . 
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Fig. 10. Time-averaged potential vorticity fields for the ex- 
periment in which only \f/ ohs is assimilated at the surface, (a) 
Layer 1 : contour interval is 5 X 10' 6 s' 1 . (b) Layer 2: contour 
interval is 2.5 X 10 -6 s -1 . (c) Layer 3: contour interval is 2.5 
X 1CT 6 s' 1 . (d) Layer 4: contour interval is 2.5 X 10' 6 s' 1 , (e) 
Layer 5: contour interval is 2.5 X 10~ 6 s' 1 . 


inside the closed contours. In this experiment, where 
the intensity of the eddy field is extremely weak, we 
establish the conditions hypothesized in the theory of 
Rhines and Young ( 1 982 ) . As predicted by that theory, 
advection of potential vorticity by the mean flow is 
able to establish, first, uniform values of q along 
streamlines. At this point the eddy flux term comes 
into play and tends to smooth the gradients between 
adjacent streamlines. Plateaus of “homogenized” q are 
thus created inside the closed contours, while the q 
gradients are expelled toward the rim of the gyres. 
Scatter diagrams computed on zonal lines connecting 
two different points along a closed streamline (Capo- 
tondi 1993) confirm in a more quantitative manner 


these considerations. They show, in fact, that in both 
layers 2 and 3, the q contours are very close to the ^ 
contours in a large part of the domain, as is the case 
for an inertial solution. 

Figure 1 3 shows the mean potential vorticity distri- 
butions for the second experiment. The upper-layer 
distribution is, again, practically unchanged with re- 
spect to the initial distribution as it was in the previous 
experiment. In layers 4 and 5, which are now in motion, 
larger deviations from the zonal contours of the plan- 
etary vorticity gradients can be observed. The second 
and third layers are not directly constrained by the 
nudging procedure and, in both experiments, they carry 
relatively significant components of the flow. Therefore, 
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Fig. 13. Time-averaged potential vorticity fields for the ex- 
periment in which the total ^ obs is assimilated at the surface, 
(a) Layer 1: contour interval is 5 X 10~ 6 $‘ l . (b) Layer 2: 
contour interval is 2.5 X 10 -6 s _l . (c) Layer 3: contour interval 
is 2.5 x 10~ 6 s' 1 , (d) Layer 4: contour interval is 2.5 X 1(T 6 
s~‘. (e) Layer 5: contour interval is 2.5 X 1CT 6 s 1 . 


The results of the assimilation experiments confirm 
this dynamical framework. The structure of the flow 
remains basically the same whether or not the eddy field 
is assimilated at the surface, thus confirming that the 
surface mean field ^ obs indeed defines the pattern of the 
mean streamlines in the subsurface layers. Consequently, 
this defines also the paths along which mean flow ad- 
vection of potential vorticity will take place. However, 
the evolution and final statistical steady state of the po- 
tential vorticity fields depend also upon the intensity of 
the eddy field. The characteristics of the mean potential 
vorticity distributions at statistical steady state are de- 
termined, in fact, by the relative strength of mean flow 
advection and eddy advection. If the eddies are very 
weak as in our first experiment, mean flow advection 


dominates during the adjustment phase. First, uniform 
values of potential vorticity are established along mean 
streamlines. The weak eddy mixing comes then into 
play by “slowly” eroding the gradients between adjacent 
streamlines and expelling them toward the rim of the 
gyre. If, on the other hand, eddy advection is comparable 
in strength with mean flow advection, potential vorticity 
will be stirred and mixed by the eddies before the es- 
tablishment of uniform values along mean streamlines 
is completely achieved. Therefore, only tongues of 
anomalous potential vorticity values can be observed in 
the final time-averaged distributions as a result of mean 
flow advection. 

The most remarkable difference introduced by the 
presence of an energetic eddy field is the more efficient 
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Fig. 14. Verification of the degree of agreement between stream- 
function and potential vorticity contours in layer 2. (a) Streamfunction 
contours: contour interval is 3000 mV. (b) Potential vorticity con- 
tours: contour interval is 2.5 X 10 6 s '. 


downward transfer of momentum. The intensity of the 
mean flow in the subsurface layers, in fact, is enhanced 
when eddies are assimilated. This is particularly evident 
in layers 4 and 5, where the mean flow can only be 
eddy driven. In fact, hardly any noticeable flow is found 
in these layers in the absence of an energetic eddy field. 

The results of this study also show the crucial role 
played by the surface mean field in shaping the global 
model circulation. Our choice of a climatological field 
as the mean component of the surface flow, which we 
have introduced in section 2 as a “working hypothesis”, 
is certainly debatable due to the long-term duration of 
the climatological dataset with respect to the duration 
of the Geosat dataset. The comparison of the assimi- 
lation results with in situ observations presented in Part 
11 confirms the limitations of our surface mean field. 
A successful assimilation procedure will thus require 
an accurate definition of the surface mean circulation. 
The mean fields calculated for the same period as the 
Geosat data in some recent studies (Glenn et al. 1991 ; 
Ezer et al. 1993) could be considered as possible alter- 
natives for the surface mean field in future studies. 
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Fig. 15. Verification of the degree of agreement between stream- 
function and potential vorticity contours in layer 3. (a) Streamfunction 
contours: contour interval is 3000 m 2 s~‘. (b) Potential vorticity con- 
tours: contour interval is 2.5 X 10~ 6 s -1 . 
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APPENDIX A 

The quasigeostrophic model formulation with N ar- 
bitrary layers is a straightforward extension of the two- 
layer case described by Holland (1978). Here we shall 
present the semidiscrete form of the equations (in 
which the vertical discretization has already been 
done). The horizontal discretization and the form of 
the finite-difference equations will not be discussed 
here. 
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The governing equations are the vorticity and in- 
terface height perturbation equations and the thermal 
wind relation: 

— VV* = J(f 4 VV,, tk) 4 -y- ( VV^-1/2 - Wjk + j /2 ) 

4 4 n: it = 1 to TV ( A 1 ) 


y4%,^)=^+n *=1,*. (A5 ) 

at 

Here q k is the sum of the relative vorticity V 2 / A ., the 
planetary vorticity / = /> 4 fiy, and the stretching term 
s = F k ' k - { (tk-\ - t£*) + /^m+iWa+i - tk)- The quan- 
tities are the inverse of the squared interfacial Rossby 
radii: 



*+ 1/2 


J(h k + 1/2. 'I'k+l/l) + Wk+\/2- 


k = 1 to N - 1 




/S 


(A6) 


Therefore, the carry information about the vertical 
(A2) background stratification. 


A* + l/2 = -A-(^ + |-^*)- < A3 ) 

£*+ 1/2 

Here integer subscripts (k) denote the vertical layers 
(A: increasing downward) in which the quasigeostrophic 
streamfunction is defined (nominally at the center of 
each of the layers), while fractional subscripts ( k 4 1 / 2 ) 
denote the interfaces between layers where vertical ve- 
locity and interface height perturbation are defined. 
The variables are the quasigeostrophic streamfunction 
(t k ) with horizontal velocity components ( u = — ty, 
v = tx), the interface height perturbation (h k + ]/2 ), 
positive upward, and the vertical velocity (wk+ 1 / 2 ), also 
positive upward. The horizontal coordinates are x 
(eastward) and y (northward), the Coriolis parameter 
is f =f 0 4 and the mean layer thicknesses are H k . 
The values of f 0 and 0 are chosen to represent typical 
midlatitude gyre values. The basic background vertical 
stratification is written in terms of the reduced gravity 
£*+ 1/2 = £Ap*+i/ 2 /po> where Ap*+i /2 is the (positive) 
density difference between layers k 4 1 and k. Fric- 
tional effects, written symbolically in (A1 ) as J k , are 
parameterized as lateral friction of the biharmonic kind 
(Holland 1978) in which 7 k = — A^tk- In addition, 
S7* includes a bottom friction, —tV 2 \p N , when k = N 
(the bottom layer). Note that the effect of the wind 
forcing equal to curl t/H u produces an Ekman 
pumping stretching tendency in the upper layer equiv- 
alent to a body force acting on the upper layer. The T k 
for k > 1 are zero. At the sea surface w l/2 = 0 and at 
a flat sea bottom w A + 1/2 = 0. The advective velocities 
at the interfaces needed in (A2) are calculated from a 
weighted average of the velocities in the layers; that is, 
tk+ 1/2 = (<**+ 1 / 2 )tk + ( I — «*+i/ 2 )^*+i» where a k + i /2 
= H k /(H k 4 H k+l ). 

The set of equations (A1 ), (A2), and (A3) express 
the conservation of the potential vorticity q k , defined 
as 

Qk = 4/4 F kJc -i(tk- 1 - M 

+ ^*,*+i(^*+i ” tk), (A4) 
and can be rewritten in the form 


APPENDIX B 

The successive correction algorithm, first introduced 
by Cressman (1959), is an iterative algorithm. Co- 
variance functions with decreasing radii are used in the 
different iterations in order to capture smaller and 
smaller scales. The procedure is similar to the one out- 
lined by Roemmich (1983) in his estimation of hy- 
drographic quantities in the Strait of Florida. The spa- 
tial covariance functions are Cressman functions (R 2 
— r 2 )/(R 2 4 r 2 ), where r is the distance between the 
data point and the analyzed point and R is the corre- 
lation distance. At each iteration corrections to the 
previous estimate are computed at each grid point ac- 
cording to the following formula: 

F°{v 4 1) 

Z »*(»)(*■?-/■?(»)) 

= FJ(») + ^ n , (Bl) 

Z Vi’v.ft') 

1=1 

where jFJ(v) is the interpolated (analyzed) value at 
position x, iteration v ; F® is an observed value at po- 
sition i ; and F f ( v ) is the estimate of the field at position 
i for the uth iteration. The weights w X i(v) involve sep- 
arate space and time factors: 

( exp[— (f* — /,) 2 /1.44r 2 ] 

+»(v) = < X (Rl - rli)l(Rl + r 2 xi ), R > r, (B2) 

1 0, r> R 

where /, is the time at observation point /, r xi is the 
spatial distance between interpolation point x and ob- 
servation point /, and T and R v are correlation time 
and space scales, respectively. Thirty-four days of al- 
timeter data centered on the time of the analysis are 
introduced with an e-folding scale, 7\ of 5 days for 
each of the analyses, performed every 2 days; T was 
chosen to be 5 days, a time much shorter than the 34- 
day search window. This value for the Gaussian e-fold- 
ing scale was intended to include altimeter data points 
from tracks west ( — 3 days) and east (43 days) of any 
particular track, while preserving ocean signals with 
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periods in the mesoscale band and longer (V. Zlotnicki 
1990, personal communication). The spatial radii R v 
are 2°, 1.5°, 1.25°, and 1° for each of the four itera- 
tions, respectively. These values were determined by 
trial and error to capture scales as small as possible 
while still producing smooth maps. 

The successive correction algorithm is not an opti- 
mal algorithm since the weights w xi in ( B 1 ) are not 
chosen so to minimize the expectation value of the 
interpolation error (Bretherton et al. 1976). The reason 
for the choice of a suboptimal algorithm is mainly as- 
sociated with computational efficiency. Due to the need 
of inverting large matrices, the optimal interpolation 
schemes involve a large computational load, especially 
when the analyzed fields are computed very frequently 
in time. Also, the assumptions underlying the deter- 
mination of the optimal weights (isotropy and ho- 
mogeneity of the statistics, knowledge of the error sta- 
tistics, etc.) are often not met in the actual implemen- 
tation of the algorithm, thus reducing its optimal 
character. 

The successive correction scheme does not auto- 
matically supply error maps for the interpolated fields. 
Interpolation error maps should be considered together 
with the interpolated fields to be able to assess the rel- 
ative reliability of different estimated values. One can 
compute error fields for the successive correction al- 
gorithm as outlined by Wunsch ( 1989). However, this 
would introduce a computational load comparable 
with the one of the optimal interpolation algorithm 
itself. Since we are not using an assimilation scheme 
that can rigorously account for the data error distri- 
bution, we delay for the moment the computation of 
the error maps. 
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